Finite-size scaling in extreme statistics 
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We study the convergence and shape correction to the limit distributions of extreme values due 
to the finite size (FS) of data sets. A renormalization method is introduced for the case of indepen- 
dent, identically distributed (iid) variables, showing that the iid universality classes are subdivided 
according to the exponent of the FS convergence, which determines the leading order FS shape cor- 
rection function as well. We find that, for the correlated systems of subcritical percolation and l/f a 
stationary (a < 1) noise, the iid shape correction compares favorably to simulations. Furthermore, 
for the strongly correlated regime (a > 1) of 1/ f a noise, the shape correction is obtained in terms 
of the limit distribution itself. 
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Extreme value statistics (EVS) has been much stud- 
ied in engineering finance and environmental sci- 
ences Q where extreme events may have disastrous con- 
sequences. Recently, interest in EVS has also been grow- 
ing in physics, e.g. in spin glasses Ufl, interface fluctu- 
ations [5j, [y] , and front propagations [7( . Unfortunately, 
the use of EVS is hampered by the cost of acquiring good 
quality statistics: EVS is derived from the extremes of 
subsets of a data set, requiring abundant data for rea- 
sonable statistics. Data analysis is further complicated 
by the fact that, while the EVS limit distribution may 
be known, the convergence with increasing sample size 
is slow. Clearly, a detailed finite-size (FS) analysis pro- 
viding the convergence rate and shape corrections to the 
limit distribution is much needed. While for iid variables 
FS studies exist in the mathematical literature [H, for 
correlated systems the convergence rate and shape cor- 
rections are known only in a few cases, such as Brownian 
motion [§]. 

In this Letter we use analytic and phenomcnological 
approaches, combined with simulations, to investigate 
FS scaling in EVS. First, we develop a renormalization 
group (RG) method, in which the limit distribution is 
a fixed point of the flow in function space of the finite- 
sample EVS distributions. Applied to iid variables, the 
approach provides an intuitive and accessible summary 
of the mathematical results for the leading FS correction, 
including the explicit forms of the shape corrections (scal- 
ing functions). Next, we consider two systems with corre- 
lated variables, namely percolation and 1/ f a signals. We 
numerically study the distribution of the largest clusters 
in subcritical percolation. While the limit distribution is 
known to be an iid problem 1(| U | , we find that even the 
FS correction fits the iid prediction well. In the case of 
the maximum statistics of 1/ f a signals, < a < 1 corre- 
sponds to the weakly correlated regime, with an iid limit 
distribution [l2|. Our simulations indicate that the FS 
properties are very close to the iid case for < a < 0.5, 
but deviations appear for 0.5 < a < 1. For a > 1, how- 
ever, the convergence becomes fast (power law) and we 



can show that, under a mild assumption, the FS shape 
correction is given in terms of the limit distribution and, 
furthermore, the order as well as the shape of the cor- 
rection strongly depends on the way the distribution is 
scaled. The paper is concluded by remarks on higher 
order FS corrections. 

The case of iid variables has been extensively studied 
fl3| . and we begin our FS study by a reinterpretation 
of the original derivation of the extreme limit distribu- 
tions 14[. Consider random variables zi, z%, ■ ■ ■ , Zn with 



parent density p(z) and integrated distribution fi(z) = 
I- p( s )ds. The maximum of the z% has the integrated 
distribution /J. N (z) and the basic observation 14| is that 
if, after an appropriate scale change z = a^x + bpf, /J. N (z) 
tends to a limit distribution M(x) as N — > oo, then the 
same limit should be reproduced for another N' = pN. 
This requirement can be cast in the form 



M{x) = [R p M](x) ee M p (a p x + b p ) 



(1) 



where the r.h.s. defines R p , which can be interpreted as 
an RG operator based on the analogy with critical phe- 
nomena. Indeed, the operation of raising to power p > 1 
and shifting by b p eliminates the irrelevant small argu- 
ment part of the parent distribution, a p rescales the rele- 
vant "degrees of freedom" , and the fixed point condition, 
Eq. {l}, determines the limit distribution. 

Eq. J]} is solved by a p = p 7 , b p = 7~ 1 (p 7 — 1) and 
M (x) = cxp [— (1 + jx) 1 ^ 7 ] , where the final scale of x and 
the position of the distribution are set by the standard- 
ization M (0) = M'(0) = 1/e. We thus have a line of fixed 
points parameterized by 7, and M(x) is the generalized 
extreme value distribution. The traditional universality 
classes are called Frechet (power decay of parent at in- 
finity), FTG (Fisher-Tippett-Gumbel, faster than power 
decay), and Weibull (power decay at a finite cutoff), and 
correspond to 7 > 0, = 0, < 0, respectively. 

In the RG picture, the FS behavior is determined by 
the action of the RG transformation on the neighborhood 
of the fixed point M(x). Thus, we consider distributions 
as M e (x) = M(x + etp(x)), assuming e is small. Stan- 
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dardization implies ^(0) = tp'(0) = 0, and the scale of e 
is set by il>"(0) = 1. Our central observation is that the 
large N behavior corresponds to the eigenvalue problem 



M e .(x) = [R p M e }{x) = MP{a p . e x + b p , e ), 



(2) 



where linearization in e is understood, a Pje , b Pi€ differ from 
the fixed point values a p ,b p determined above to 0(e), 
and the eigenvalue is A = e'/e. Using the fixed point 
relation (JTJ), we obtain [\/ a p )ip" {x) = -0" (a p x+b p ) whose 
solution with p- independent ip(x) reads as 



ip(x) 

\ = p 1 



(1+7x)t'/t +1 - ( 7 '+ 7 )x-l] / 7 '(7'+7) (3) 

(4) 



Thus we see that, for a given universality class parame- 
terized by 7, a new parameter 7' emerges characterizing 
the eigenvalues and eigenfunctions of Eq. ([2]). Note that 
shape corrections equivalent to ip have been obtained by 
direct methods in the mathematical literature [H, [H| . 

In order to link the RG result with the N dependence, 
we write e = ejv and use ^ to find e' = e p N = p 1 ej\r. 
Assuming a power form one obtains 



e N oc N 1 



(5) 



or, more precisely, 757%^ ■ ^ 7'- Thus 7' is the FS con- 
vergence rate (stability implies 7' < 0). To find en and 
thus 7' for a given parent, fx(z), we study the integrated 
distribution function p, N (z) with the shift and scale pa- 
rameters &at — /i(lniV), ajy = /i'(ln AT) expressed through 
h(y) = )jr 1 {e~ e ). Close to the fixed point one has 

M N (x) = n N (a N x + b N ) w M(a; + e N ip(x)), (6) 

and differentiating — ln[— lnMjv(a;)] twice at a; = gives, 
to leading order, da^ /db]\r — > 7, and at next order 



£JV = 7 — da,N /db^ ~ A^ 7 



(7) 



The convergence rate 7' is now determined and the per- 
turbation function t/>(cc) follows from (jSJ. This gives a 
practical meaning to the results from RG theory. 

For data analysis it is convenient to represent the FS 
correction with zero mean (x) and unit variance ex 2 (finite 
for 7 < 1/2) by using the variable y = (x— (x))/cr x . Here 
we consider the FTG class (7 = 0) with limit distribution 
M(°\y) = e~ e ' v+ ) where a = it / \f&, b = 7^, the latter 
being Euler's constant. Writing Mn(jj) ~ M^(y) + 
ejvM^^y), we have for the correction 

M^(y)= P {0 \y)[e~<' (a y + V + ay + /3] / a-/' 2 , (8) 

where P (0) (y) = M^{y)', a = T(l - j') b +*^-^ , with 
#(z) = T'{z)/T{z) and = -T(l - 7'). For 7' = 0, 
Eq.© becomes [(,(z) denotes Riemann's zeta function] 




M^(y)= P {0) (y)[a 3 (y 2 



2C(3)y]/2a 2 . (9) 
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FIG. 1: Finite-size analysis for the largest clusters in subcrit- 
ical percolation. The prediction from iid theory (solid line) is 

$i(y) = M m (yY f rom E q- ©• The simulation results (dot- 
ted line) were obtained with occupation probability p — 0.25 
and system sizes L — 500, 1000. 

We illustrate the above results on FTG class parents 
with the commonly found asymptote 1 — fi(z) oc e~ z / z 9 
{5 > 0). Using Eq. (7}, for 5 ^ 1 we have e N w 
(<5 - l)/(8hiN), so i = 0. For 5 = 1, 9 = (expo- 
nential distribution) one finds cn ~ 1/2 AT and 7' = — 1 
while, for 5 = 1 and 6* 7^ 0, we have ew ~ — 0/ln 2 A/', so 
again 7' = 0. Thus generically 7' = 0, with FS shape 
correction given by © and the perturbation decaying 
logarithmically. Faster, 1/AT, convergence is seen only 
for 5 = 1, 9 = where Eq. fl8|) applies with 7' = —1. 

As an application, we studied the FS corrections to 
the distribution of the largest cluster size on a square 
lattice in subcritical site percolation, p < p c , where 
p c 0.592 ... is the critical occupation probability. Due 
to the finite correlation length, clusters in a large system 
are nearly independent and the size distribution of the 
largest obeys FTG, provided the inherent discreteness of 
the problem is treated appropriately T3, 11|. It remains 
an open question, however, whether the FS corrections 
can also be described by iid theory. To answer this ques- 
tion, we first note that the asymptote of the distribution 
of the cluster size s is s _1 exp(— s/s^) [l6| . where is the 
cut-off size. This asymptote corresponds to 5 = 9 = 1 
in the example of the previous paragraph, so © gives 
the FS scaling function and sw — 1/ln 2 N. To com- 
pare theory with simulation, we collected statistics for 
the largest cluster in systems of sizes L = 500 and 1000 
in an ensemble of sa 10 7 runs, resulting in a relatively 
smooth histogram. The shape correction was then ob- 
tained by subtracting the empirical histogram from the 
FTG distribution, multiplied by In 2 N, where AT is the 
average number of clusters. The result is compared with 
the iid theory in Fig. [TJ As can be seen the curves match 
surprisingly well, suggesting that the iid theory is also 
relevant for the FS corrections. 

Next, we treat correlated time signals, h(t), and study 
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FIG. 2: FS shape correction to the maximum distribution 
in l/f a noise with a = 0.4. The solid line is the iid theory 
$i(y) = ~M {1) (y)'/2 from Eq. flSJ, while the others are simu- 
lations for system sizes TV = 2 5 , 2 7 , 2 9 , 2 1 , 2 15 , with minimum 
decreasing in this order. The inset shows the skewness of the 
maximum distribution for various sizes slowly converging to 
the FTG value «;§ w 1.14. It demonstrates that FS properties 
may qualitatively change near a « 0.5. 

the FS effects on the distribution of their maxima, h m . 
To have control on correlations, 1/ f a noise is chosen, 
where the Fourier amplitudes are independent Gaussian 
variables, with variance f~ a , and uniform, random phase 
171 ]. For a = 0, the process is white noise, for < a < 1 
it is stationary with correlations decaying as i Q_1 , while 
for a > 1 the fluctuations of the signal diverge as t — > oo. 
The distribution of h m has recently been studied inten- 



sively, see [18| and references therein. The main features 
to be recalled here are that the FTG distribution applies 
for < a < 1 [l2| while, for a > 1, nontrivial distribu- 
tions emerge whose shape depends on boundary condi- 
tions and the reference point from where the maximum 
is measured. Here we concentrate on the FS correction of 
the distribution of maxima, and to be specific, the max- 
imum is measured from the mean of a periodic signal. 
First, consider a < 1. Fig. [5] shows the FS shape cor- 
rection for a = 0.4 together with the iid prediction. The 
latter comes from a Gaussian parent, thus ejv ~ l/21niV, 
i = 0. The theoretical curve is $i(y) = -M^(y)'/2 
by Eq. ([9]) , corresponding to a In N magnification factor 
in the simulation curves which visibly approach $i(x). 
A similar approach can be seen for all a < 0.5. The in- 
set, showing the skewness for finite systems, also suggests 
that the leading FS correction may be described by the 
iid theory for a < 0.5. This conclusion goes beyond what 
we experienced in percolation: there correlations had a 
finite cutoff, while here correlations decay like a power. 
We now turn to a > 1, where (h m ) diverges as (h m ) ~ 
and the approach to the limit distribution 
improves from logarithmic to power-law. This effect can 
be seen in our simulations as well as in the exact results 
jtj l20| for random walks (a = 2). The limit distribution 
for a — 2 is given by the Airy distribution &Ai(z) with 



z = h m /\f~N while the first correction to scaling is 

*{z)k*m(?)+* , m (z)I>/2N . (10) 

The simplicity of the above result calls for simple expla- 
nation. Indeed, Eq. (|10p follows from the assumption 
that the shape of the distribution relaxes faster than its 
position. To see this for arbitrary a > 1, we note that 
the n-th cumulant of h m scales for large N 
with 8 = (a — l)/2 [Tij . Next, we write the corrections 
to scaling of 



(11) 



and assume that cj„ > lui for n > 1. This assumption 
implies that the shape of the function relaxes faster than 
its position. Introducing now the scaled variable z = 
h m /N e and expanding the cumulant generating function 
of h m in 7V~ Wl yields the scaled distribution function 
$ N (z) = N e P N (N e z) to first order as 



$ N (z) w <f>(z) - k\<S>'{z)N- 



(12) 



where $(z) is the scaling function in the N — > oo limit. 
For a = 2 one has cJi = 1/2 and k\ = —l/\/2 thus 
Eq. (fit)]) is recovered. 
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FIG. 3: Demonstration of the increasing speed of convergence 
to the limit distribution, $, for a = 2. Results for system 
sizes N — 16, ... , 16384 are shown using scaling variables x = 
hm/{h m ) and y — (h m — (h m )) / \f~Ka in the upper and lower 
panels, respectively. $i is the shape correction function. 
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The choice of scaling variable may change both the 
rate and amplitude of the correction term. E.g. a natural 
choice is to scale by the average [x — h m /(h m ), <&n(x) = 
(h m )PN((h m )x)]. It yields the same rate of convergence 
but it does alter the scaling function of the first order 
correction. Indeed, using x = h rn / (h m ) and expanding 
the cumulant generating function in l/N^ 1 results in 

$ N {x) « - («!/«?) [(x - l)$(x)]' N~ Ul . (13) 

The limiting function for a = 2 and its correction 

$i(x) with an amplitude 0.1 are shown in the upper panel 
of Fig. [3] Since $i(x) has two zeros, &n(x) is nearly 
fixed at two points, so the main correction comes from 
the central weight shifted to the tails of the distribution. 

If the main FS effect is due to (h m ) then the use of the 
scaling variable y = (h rn — {h m ))/y/K2 eliminates these 
corrections and, as seen in Fig. [31 improves convergence 
dramatically. The shape correction may be calculated by 
assuming that u) n > uji for all n > 2. This means that the 
leading FS correction comes from K2 and the cumulant 
generating function to first order in 1 /N U2 yields 

$jv(y) « #G0 - (k!/2kS) [y$(y) + &(y)]' . (14) 

As seen, the shape correction can again be expressed 
in terms of the limit distribution. The scaling function 
$i(y) displayed in Fig. [3] has three zeros which restrict 
possible deviations from the limit distribution to higher 
order. In addition, we found numerically o>2 ~ 3/2, re- 
sulting in such a fast convergence that the curves with 
various N-s cannot be distinguished within linewidth. 

We can make an intuitive proposition for the u> n -8, 
based on an analogy with the cumulants n n {w2) of the 
roughness u>2 of l/f a signals [T3|- There, the large N 
asymptote can be rewritten as K n (w2) ~ ( w 2) n (1 — 
b n (W2) n °^^ a ) ■ Now assuming the same exponent 
(1 — not)/ (a — 1) for K n (h m ), and reverting to the N 
dependence we get Lo n (a) = (na — l)/2 for the FS ex- 
ponent. For Brownian motion we recover 0*2(2) = 3/2, 
in accordance with simulations. Furthermore, the crite- 
rion Lo n < LOn+i is satisfied, so u>2 (a) = a — 1/2 is the 
candidate for the FS exponent. Since u) n {a) increases 
with a, the convergence is expected to improve for larger 
a. Indeed, our simulations for a = 4 show that the same 
convergence as in the lower panel of Fig.[3]can already be 
obtained in the x = h m / (h m ) scaling since cJi(4) = ^2(2). 

So far we have considered the leading FS correction. 
It is natural to ask about higher orders, especially in 
the FTG class with typically slow, logarithmic conver- 
gence. Higher order calculations are possible, which we 
just illustrate here by an arbitrary order result for the 
parent distribution fx(z) = 1 — e~ z , which, for S = 2, is 
the Rayleigh distribution, the basic proposition for the 
statistics of wave crests in ocean engineering |2l|. For 



appropriate choice of x — (z — 6jv)/ffljv as 



Mn(x) ?s exp |- 



exp 



-xH 



In TV, 



)]} + O(l/A0,(15) 



with H(u) = [(1 + u/S) —l]/u. Remarkably, all logarith- 
mic orders sum up to a scaling function in the variable 
xj In N. For a general parent distribution, the second or- 
der calculation has been carried out by a direct method 
[lit . Inspired by that, we have worked out an algorithm 
to arbitrary orders which will be presented elsewhere. 

Finally, we note that the RG approach can potentially 
be extended to the study of EVS in correlated systems. 
In cases where the limit distribution is known, the FS 
corrections may be clarified, while in less explored sys- 
tems it may help in finding the limit distribution itself. 

This work was supported by the Hungarian Academy 
of Sciences (Grant No. OTKA K68109). NRM acknowl- 
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<5 = 2, the distribution of the maxima is obtained with 
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